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We present a hybrid lattice Boltzmann algorithm for the simulation of flow glass-forming fluids, characterized 
by slow structural relaxation, at the level of the Navier-Stokes equation. The fluid is described in terms 
of a nonlinear integral constitutive equation, relating the stress tensor locally to the history of flow. As an 
application, we present results for an integral nonlinear Maxwell model that combines the effects of (linear) 
viscoelasticity and (nonlinear) shear thinning. We discuss the transient dynamics of velocities, shear stresses, 
and normal stress differences in planar pressure-driven channel flow, after switching on (startup) and off 
(cessation) of the driving pressure. This transient dynamics depends nontrivially on the channel width due 
to an interplay between hydrodynamic momentum diffusion and slow structural relaxation. 
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I. INTRODUCTION 

The flow of glass forming fluids is characterized by 
an interplay of slow, collective structural relaxation and 
flow-induced relaxation. In many applications in particu¬ 
lar in colloidal suspensions, the structural relaxation rate 
1/r is much larger than the flow rate 7; hence the relevant 
Peclet number Pe = yr 3> 1. This leads to pronounced 
nonlinear-response effects such as shear-thinning: the ef¬ 
fective viscosity of the fluid decreases strongly with in¬ 
creasing 7. This regime opens close to the glass tran¬ 
sition, where r —>■ 00, and hence Pe 3> 1 even if the 
shear rate is still slow compared to the individual-particle 
short-time relaxaton rate I/tq, such that the bare Peclet 
number Pep =_7'n3 ^ 1- This is the regime of nonlinear 
glassy rheologyiJ. 

In transient dynamics, viscoelastic and shear-thinning 
effects combine: on time scales t :S> tq but t <C I/7, the 
system effectiely behaves as a transiently frozen amor¬ 
phous structure characterized by some elastic modulus 
Goo • This was first recognized by MaxwelP for the linear- 
response regime (Pe ^ 1). The simplest fluid model 
describing this phenomenon, called Maxwell model in at¬ 
tribution to him, can be written as 

<^xyit)= dt', (1) 

Jo 

where simple-shear flow 7 = dyV^ is assumed to start at 
t = 0 from a stress-free equilibrated state, and cr^yit) is 
the shear-stress component of the Cauchy stress tensor 
(T(t). For t T, the model gives the stress-strain re¬ 
lation of an elastic Hookean solid, a{t) ^ Goolto where 
7tt' = /// 7(5) ds is the accumulated shear strain. For 
t T and constant shear rate, a{t) ^ 777 recovers vis¬ 
cous Newtonian flow with a shear viscosity rj = GooT 
given by the so-called Maxwell relation. 
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Equations such as Eq. ([^ are constitutive equations for 
continuum mechanics: the macroscopic flow field v and 
its gradients k = (Vv)^ are determined by the Navier- 
Stokes equatiorP, 


dtgv -j- V • {gw) = gf - VpV ■ er , (2) 


where / is a given external force density, and p is the 
thermodynamic pressure. The mass-density held g obeys 
a continuity equation, dtg + 'V ■ gv = 0. We assume in 
the following incompressible how, V • v = 0. The stress 
tensor a expresses microscopic friction effects that need 
to be described by a constitutive equation that expresses 
cr in terms of the how helds again. 

Recent theoretical work based on a formalism called in¬ 
tegration through transients (ITT), developed by Fuchs 
and Cates in the context of driven colloidal huidd^, gives 
a framework to derive constitutive equations for non¬ 
linear glassy rheology from microscopic theory. Assum¬ 
ing how to remain homogeneous at least on mesoscopic 
scales, one gets, schematically, 

o-struc(0 = f [-dt'B{t,t')]G{t,t',[K])dt' , (3) 

Jo 

where the subscript “struc” recalls that this is only the 
structural-relaxation contribution to the stresses (akin to 
the purely polymeric stress contribution in polymer rhe¬ 
ology). G{t,t', [k]) is a nonlinear dynamical shear mod¬ 
ulus that depends on the whole history of deformation 
gradients in a suitable way that honors the invariance 
of stresses under rigid-body motions (called material- 
frame indifference in continuum mechanics). Replac¬ 
ing it with a how-independent exponential recovers the 
Maxwell model. B{t,f) = E{t,t')-E^{t,t') is the Finger 
tensor, related to the deformation tensor 


E{t, t') = exp^ 


k{s) ds 


( 4 ) 


Here exp^ denotes a time-ordered exponential where all 
products of K are sorted such that earlier times appear 
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to the right. The appearance of in Eq. ^ 

generalizes the scalar shear rate in Eq. Q to arbitrary 
incompressible flow geometries. 

A natural feature of constitutive equations such as 
Eq. (§ is their integral nature. Other than Eq. 0 , it 
cannot in general be reduced to a differential equation in¬ 
volving only time-local differential operators. This is in 
particular relevant close to the kinetic glass transition, 
where the dominant feature in the dynamics is the slow 
structural relaxation that provides long-lasting memory 
effects. Within the mode-coupling theory of the glass 
transition (MCT), the dynamical shear modulus is deter¬ 
mined from the solutions of a set of integro-differential 
equations for density correlation functions that contain 
the effects of flow advectiorP. As a key feature, the his¬ 
tory integral does not have a predefined natural cutoff, 
but rather it extends arbitrarily far back in time, in a 
way that depends on the state point and the flow his¬ 
tory. This in essence reflects a lack of separation of time 
scales into short microscopic ones and potentially slow 
hydrodynamic ones. 

To study the flow of glass-forming fluids, it is hence de¬ 
sirable to develop continuum-mechanics solvers that can 
combine the Navier-Stokes equation with integral-type 
constitutive equations such as Eq. (^. In this paper, we 
describe such a sc heme, based on the lattice Boltzmann 
(LB) algorithirP^ to solve the Navier-Stokes equation in 
the low-Mach number limit. 

The LB algorithm is a kinetic scheme based on lattice- 
node densities that evolve according to collision-and- 
streaming rules taylored to reproduce in the continuum 
limit of vanishing lattice spacing and time-step length 
the Navier-Stokes equation for a Newtonian fluid. We 
base our work on our recent LB algorithm incorporat¬ 
ing non-Newtonian stresses through a modified collision 
rul^. This algorithm is here extended to a full “hybrid- 
LB” scheme combining LB steps for the Navier-Stokes 
equation with an integral-equation solver keping track of 
the full flow history in Euler coordinates. Previous ap¬ 
proaches extending LB by additional lattice populat ions 
or forcing terms have focused on linear viscoelasticitjl^Ufll 
and nonlinear models that can be expressed in terms 
of differential constitutive equationJSHi^ for a review 
see also Ref. m Also, modified collision rules can be 
used to simulate the nonlinear rheology of emulsions, ex¬ 
ploring the flexibility of the LB algorithm as a kinetic 
schem^i^IIIl 

Earlier hybrid-LB schemes have coupled the 
LB algorithm with finite-difference solvers for differential 
constitutive equations, usually entering non-Newtonian 
stresses in terms of a body forc^^oHlIl Qur hybrid-LB 
scheme is new in its focus on integral constitutive equa¬ 
tions with large relaxation times. For the treatment 
of integral constitutive equations combined with finite- 
element and finite-volume algorithms, see Refs. [23] and 
123 for a recent LB-based finite-volume formulation, see 

also Ref. |23 

As a specific application, we implement a nonlinear 
generalization of the Maxwell model that combines the 


effects of shear-thinning and viscoelasticity. We use this 
model to study generic transient effects of the startup 
and cessation of pressure-driven flow in a planer (2D) 
channel. 

The paper is organized as follows: we introduce in 
Sec.[T^the integral nonlinear Maxwell model. In Sec. |III[ 
we briefly describe the LB algorithm based on Ref. | 8 | 
and our integral solver. Section presents results for 
pressure-driven channel flow, followed by a concluding 
Sec. El 


II. NONLINEAR MAXWELL MODEL 

T he key interplay of mechanisms expressed by the ITT- 
MCtREEZIMI jg tJiat of slow structural relaxation of rate 
T~^, and of flow-induced relaxation of rate | 7 |/ 7 c (where 
we introduce 7 c as a model parameter controlling the ef¬ 
fectiveness of strain in breaking nearest-neighbor cages). 
If Pe 1 » Peo, the integral in Eq. 0 is cut off times 
t' <t — 7 c/| 7 |, leading to a decrease of the effective vis¬ 
cosity, i.e., shear thinning. 

In the Maxwell model, this can be incorporated by 
letting (in steady state) >->• + | 7 |/ 7 c- For general 

time-dependent incompressible flow, a plausible choice is 
to set, in Eq. 

Git, a, M) = ^ (5) 

For the flow rate, we set 7 ^ = //d = (1/2) trZ)^ (for 
incompressible flow), where D = k, -V is the sym¬ 
metric velocity-gradient tensor. The quantity 7 c is set 
to 1/10 in our numerical calculations; this follows previ¬ 
ous ITT-MCT investigations of colloidal glass former^^ 
where it was related to the typical fraction of a parti¬ 
cle diameter each particle can be sheared before cages 
break (agreeing with an empirical criterion for the melt¬ 
ing of solids given by Lindemann). Our integral nonlin¬ 
ear Maxwell model (nlM) reproduces qualitative features 
found in experiment and ITT-MCT for time-dependent 
nonlinear glassy rheology, e.g., for large-a mplitu de oscil¬ 
latory shear or creep under imposed stresd^^*^. It also 
describes the discontinuous emergence of a finite yield 
stress at the glass transition, t ^ 00. It should however 
be kept in mind that Eq. <0 represents a gross oversim¬ 
plification of the ITT-MCT dynamics. The identification 
of II as the local shear rate is ad-hoc. It is the sim¬ 
plest material-frame indifferent choice that is linear and 
respects the symmetry of flow reversal (under which the 
relaxation rate must remain positive). A notable feature 
of our model that is faithful to the microscopic ITT-MCT 
constitutive equation, is that the dynamical shear mod¬ 
ulus defined by Eq. 0, breaks time-translational invari¬ 
ance for non-stationary flow. For this reason it cannot 
be reduced to a differential constitutive equation. 

Equation 0 together with Eq. 0 reduce to the well- 
known upper-convected Maxwell model (UCM)P^ in the 
case where Pe <C I (such that the second exponential in 
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Eq. (§ can be approximated by unity). Although the 
UCM can be written in differential form, it will serve 
as a useful test of our algorithm since analytical solu¬ 
tions are available for the transient dynamics (see Ap¬ 
pendix]^. A different generalization of the UCM to non¬ 
linear rheology in terms of a differential equation is the 
White-Metzner modeP^ it corresponds to replacing the 
shear-rate-dependent expression in the second exponen¬ 
tial of Eq. eqrefnlmax by the accumulated strain. Which 
is closer to true ITT-MCT may actually depend on the 
type of flow considered. ITT-MCT keeps a much more 
complicated strain dependence that is not easily cast into 
a simple form for the shear modulus. 

The Maxwell model only addresses structural relax¬ 
ation on times large compared to tq. On this short-time 
scale, the dynamical shear modulus decays from its in¬ 
stantaneous value to the Maxwell plateau modulus Gao- 
We are not concerned with this regime here, and simply 
assume this process to provide a Newtonian background 
viscosity rjoa = Goo To that is shear-rate independent. For 
colloidal suspensions, this may be thought of as a crude 
model of the solvent viscosity, ignoring the flow-induced 
hydrodynamic interaction effects. We hence set 

cr(t) = cr struc(t) + r]ooD{t), (6) 

III. LATTICE BOLTZMANN METHOD 

The lattice Boltzmann method is a fast and versatile 
tool to solve the Navier-Stokes equations. For a review, 
we refer to Ref. [71 Considering a regular rectangular spa¬ 
tial grid with lattice spacing Sx, the LB scheme evolves 
a set of lattice-density distributions rii are associated to 
a finite number of velocities that represent streaming 
from a node to, usually, the nearest and next-to-nearest 
neighbors. The lattice distributions are evolved over a 
time step St by a collision-and-streaming rule 

ni{r CiSt, t + 6t) = n* (r, t) 

= ni{r,t) + A-t[n{r,t)]+Fi. (7) 


where the equilibrium distributions are given by 

u )■ <"> 

The lattice weights a° = 4/9, = 1/9, = 1/36, and 

the speed of sound Cg = c/Vi are chosen to reproduce 
the flow of a Newtonian fluid if the forcing term Fi is 
set to zero. In this case, the Newtonian shear viscosity 
is given by ? 7 n = {St)pcVTLB — 1/2). Note that the equi¬ 
librium distribution depends only on the fluid density p 
and velocity u, but not on the flow rate. 

To model non-Newtonian stresses cr"^, we setP 


F = a‘^ 


-1 

2c4rLB 


_flJN / Z JT \ 


fex 

Ja ^10 






2ci 


( - {St){dtSp)UaUj3 + {upfV + UafV)) | • (10) 


As the non-Newtonian stresses are not in general trace¬ 
less, we split into its traceless part denoted by 
an overbar, and a non-Newtonian pressure contribution 
jpiiN _ _i/ 2 cr/;J^, where cr"^ denotes the trace of the 
non-Newtonian stress tensor. In incompressible flow, the 
equation of the state of the D2Q9 model, po = pc/, re¬ 
lates this extra pressure term to a small variation in the 
density which is implemented via 


dtSp{t) = • (11) 


The hydrodynamic density and momentum fields are re¬ 
covered from the relation^ 

Pir,t) = Y^n^ + jdtp, (12a) 

i i 

pu{r,t) = ^CinV = + . ( 12 b) 


The collision operator implements relaxation towards 
a set of equilibrium distributions nV that are chosen such 
that for a specific lattice, the desired continuum limit 
emerges as 5x, St 0. The term Fi is used to model 
external forces or, in our case, a non-Newtonian part of 
the stress tensor. 

We consider a two-dimensional grid for simplicity; the 
extension to 3D is straightforward. The velocity set is 
chosen according to the standard D2Q9 model incor¬ 
porating nine lattice velocities: Cq = ( 0 , 0 ), Ci. ,4 = 
(± 1 , 0 )c, ( 0 , ±l)c, and c^...Vt = (±l,±l)c, in units of 
the lattice velocity c = Sx/St. For the collision operator, 
a single-relaxation time BGK model is employed, 

Ai = -—{ni-nV) ( 8 ) 


For the justification of this LB scheme, we refer to Ref. HI 
where a standard Chapman-Enskog expansion was used 
to demonstrate that the continuum limit of our scheme 
is indeed the Navier-Stokes equation supplemented with 
a non-Newtonian stress contribution. 

For the application to planar channel flow considered 
below, we assume the flow to remain translational in¬ 
variant for computational efhcency. Generalized periodic 
boundary condition^^ are used to fix a constant pressure 
step between the inlet and outlet of the periodic lattice. 
Under these conditions, only the terms quadratic in the 
velocities need to be kept in Eq. (10). The definitions 
of the hydrodynamic fields, Eq. (12), then reduce to the 
standard ones discussed in the LB literatur^. 

To implement constitutive equations of the form of 
Eq. (§, we keep track of the Finger tensor B during the 




















4 


simulation, introducing a time-history grid for the inte¬ 
gration of the constitutive equation at each LB lattice 
node. Since in the regime of interest for nonlinear glassy 
rheology, the structural relaxation time t can be orders 
of magnitude larger than tq, the time integral in Eq. ^ 
extends backwards over a potentially large time span. 
To deal with this, we employ a quasi-logarithmic mem¬ 
ory layout consisting of B blocks (labeled b = 1,... B) of 
equidistant lattices with C grid points of fixed time step 
Stf,. The time steps are doubled from block to block, 
as time extends backwards from t io t' < t in the in¬ 
tegration: = 2Sti,, identifying the smallest step 

size as that of the LB solver, Sts = St. This quasi- 
logarithmic grid assumes that the function G{t,t',[K]) 
entering Eq. (§ varies slowly for large t — f; this is indeed 
a feature of both our nlM model and the full ITT-MCT 
if the time-dependence of the flow rate itself is not fast. 

Instead of computing the time derivative of the Finger 
tensor directly, we save the velocity gradient tensor k, and 
the discrete contributions exp{Kdt) to the deformation 
tensor 


_ ^K{t—6t)5t ^n(t — {C—l)St)St^K,{t — CSt)St 

...e'^^t')StB-x ^ (13) 

As the LB scheme steps forward in time and a time- 
integration block is filled, the two oldest entries are mul¬ 
tiplied and moved to the underlying block with a time 
step twice as big. Within lattice accuracy, this proce¬ 
dure keeps the exact value of the deformation tensor. 
For the velocity-gradient tensor k a further approxima¬ 
tion is needed: we keep the averaged tensor over the two 
oldest entries when transferring them to the next block 
backwards. The time derivative of the Finger tensor is 
then evaluated according to 




(14) 


At each LB lattice node, the integration of Eq. @ 
can then be performed by a suitable integration scheme. 
For the slowly varying functions we expect on physi¬ 
cal grounds, a simple trapezoidal rule is sufficient. Our 
hybrid-LB scheme is particularly adapted to these situ¬ 
ations where the constitutive equation is given in terms 
of Euler coordinates (as exemplified by the appearance 
of the Finger tensor, instead of the Cauchy-Green ten¬ 
sor). Thus, we do not need to keep track of the flow- 
advected movement of Lagrangian material points; at the 
expense of needing to evaluate the time-ordered exponen¬ 
tial Eq. (13) based on the generator k of the nonlinear 


deformation. 

In the following we will test the scheme for a simple 2D 
setup that is translation-invariant in the flow direction. 
Assuming the flow to remain laminar, we considerably 
increase the efficiency of the computationally demanding 
algorithm (approximately 10® lattice-node updates per 


second using 1 core of an Intel Core i5-3470S CPU) by 
evaluating the non-Newtonian stress only in the central 
column of the lattice and relaying this extra stress along 
the symmetry axis. We have checked in separate simula¬ 
tions that this does not affect the results. Unless stated 
otherwise, we use tlb = 0.9 and a grid of 100 x 20 lat¬ 
tice nodes (depending on the Maxwell relaxation time r). 
For the integration of the flow history, we have chosen a 
block size of C = 128. The velocity gradient tensor k is 
evaluated using a second-order finite difference scheme. 
The algorithm is implemented in the open-source lattice 
Boltzmann code Palabo^^. 

For the simulations, we choose parameters as follows: 
a pressure drop along a channel of length L is consid¬ 
ered that is comparable to the Maxwell elastic modulus, 
Ap/L = Gao', this is typical for soft-matter fluids, where 
Gao = 0(1 Pa). Our model is that of a yield-stress fluid, 
and this pressure difference is sufficient to flow-melt the 
glassy state modeled by r = oo. Times are measured in 
units of the microscopic relaxation time tq (on the order 
of 1ms for typical colloidal fluids); the relevant param¬ 
eter expressing the viscoelasticity of the system is then 
the relative slowdown of structural relaxation, 9 — t/tq. 
Unless stated otherwise, calculations are performed for 
9 = 10. The dynamics does not change qualitatively for 
higher values of 9. 


IV. RESULTS 
A. Steady State 

We briefly discuss the stationary flow profiles. For 
translational-invariant flow, the model defined by 
Eqs. ([^ and ^ can be solved analytically to give 

o-ss = f?oo (« + , (15a) 

n>l 


with 



Rereading this as stress-strain-rate relation at arbitrary 
time t defines the instantaneous nonlinear Maxwell model 
used in Ref. |H1 Such instantaneous constitutive equations 
(that relate the stresses cr(t) to the strain rates K{t) at 
the same time) are numerically much less demanding to 
implement. However, they only account for steady-state 
shear thinning and not for viscoelasticity. To check that 
the integration scheme used here for the flow-history inte¬ 
gral converges to the correct steady state, we compare in 
Fig.g the velocity, shear-stress, and normal-stress pro¬ 
files of the two models with the analytical solution for 
pressure-driven channel flow given in Ref. [H Both LB 
results match the analytical solution perfectly for the 
value of r = IOtq. This value is chosen as a moderately 
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FIG. 1. Profiles of the steady-state velocity (boxes; left y- 
axis), shear stress (triangles), normal-stress difference (circles, 
multiplied by 20), and pressure (diamonds, multiplied by 20; 
right y-axis) for flow in a 2D channel of width 2H (driven by a 
pressure gradient Ap/ {2H) = Goa), for the nonlinear Maxwell 
model, Eq. with r/ro = 10. Open symbols (left half) 
are for an instantaneous model, Ref. |51 closed symbols (right 
half) for the present model, obtained from our LB algorithm; 
lines are analytical solutions. The Newtonian (dashed line) 
and glassy (dotted line) profiles are shown for comparison. 
Arrows indicate the positions used in later plots. 


viscoelastic case (keeping the numerical effort in solving 
the history integrals moderate), that already represents 
many features of the glassy solution (r —oo, shown as 
dotted lines in Fig.[^. 

Figure demonstrates typical effect for the channel 
flow of viscoelastic shear-thinning fluids: for pressure gra¬ 
dients per unit length comparable to the shear modulus, 
it is advantageous for the fluid to form high-shear re¬ 
gions near the walls and a co-moving low-shear “plug” 
in the center. This is explained by the finite yield stress 
that arises at the glass transition. In the center of the 
channel, the shear stress remains below the yield stress, 
so that no homogeneous flow is possible there. For our 
model, the yield stress has a shear-stress component of 
o'y,a;y = Goo7c- The tensorial structure of the model 
includes normal stresses (circles and diamonds in the 
figure). They are quadratic in the shear-rate due to 
material-frame indifference of the stress tensor. In the 
flowing region outside the plug, the model predicts a 
positive normal stress difference, ~ ’^yy > 0- This 
causes a force towards the channel center, which in the 
incompressible fluid is balanced by an increase in the lo¬ 
cal density. In colloidal suspensions, this can be a driving 
mechanism for particle migration to the centei^^. 


B. Transient dynamics 

We now consider the transient dynamics of the channel 
flow upon applying and removing a sudden driving pres¬ 
sure step. Figure [^presents an overview of the transient 
evolution for both cases, for the velocity (top), shear- 
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FIG. 2. Velocity u (top), shear stress a^y (middle), and 
normal-stress difference a^x — Oyy (bottom) evolution after 
startup (left) and cessation (right) of pressure-driven channel 
flow of a nonlinear Maxwell fluid (relaxation time 6 = 10), 
evaluated at position di /2 = v/H = 0.51 in the channel of 
width 2H, and normalized by the respective steady-state val¬ 
ues. Curves are shown for different channel widths corre¬ 
sponding to twaii/To = 0.01, 0.1, 0.5, 1, 2, 4, and 8 (heavy 
lines); light grey lines mark values of twaii equidistantly spaced 
in between. 


stress (middle), and normal-stress difference (bottom) at 
a quarter-width of the channel (position dif 2 marked in 
Fig. a Curves are shown for several channel widths. 
While for wide channels, the velocity u{t) increases mono- 
tonically after application of the pressure gradient, the 
evolution to the steady state is nonmonotonic for narrow 
channels. 

In a Newtonian fluid, the channel diameter defines the 
only relevant time scale of the transient dynamics. As 
long as the flow is purely laminar, transverse-momentum 
diffusion across the channel sets the time scale over which 
information of the boundary conditions is transmitted to 
the center. This time scale is inversely proportional to 
the Newtonian viscosity rjao = Goo 7), and read^^ 


^wall 


4H^p 

TT^GooTo ' 


( 16 ) 


As seen in Fig. twaii ~ 1 separates the regime of nar¬ 
row channels from the regime of wide channels for the 
evolution of the startup velocity. This can be retional- 
ized as follows: At short times after startup, t tq, 
the fluid still behaves as a Newtonian fluid with vis¬ 
cosity ryooi since the structural-relaxation contribution 
to the stresses given by Eq. ^ is still negligible. For 
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twaii ^ To, this time is sufficient to transiently build up 
the parabolic velocity profile of the Newtonian fluid in 
the channel. Only at times t twain the flattened “plug¬ 
like” velocity profile of the non-Newtonian fluid will be 
established. Hence, the transient velocity first increases 
towards the larger Newtonian steady-state value, before 
it decreases again to the lower non-Newtonian one. If 
twaii To, the initial high-shear Newtonian profile is not 
established, and the velocity monotonically increases to 
the non-Newtonian steady state. 


For the evolution of the shear stress (Txyit) (middle 
left panel of Fig. [^, twaii marks the crossover between 
the initial rise and a slower approach to the steady state. 
The normal-stress differences shown in the bottom panel 
of the figure display a more complex evolution towards 
the steady state, with overshoots discernible for twaii ^ 
To, and a monotonic increase for the narrowest channels 
shown. This pattern depends on the position y/H across 
the channel, as will be discussed below. 

Starting and stopping flows are symmetric for New¬ 
tonian fluids, in the sense that u{t) after startup and 
u{t) = Uss — u{t) after cessation are identical (where Uss 
is the steady-state value(P^. The same symmetry holds 
for linear viscoelastic models such as the UCIVPS. It is 
broken for nonlinear constitutive equations, since there 
the stress-strain-rate relations after cessation (with no 
flow present) differ from that after startup (with flow 
present). This can be seen in the right panels of Fig. 
the velocity after cessation decays nonmonotonically to¬ 
wards zero even for twaii ^ 1- The asymmetry is par¬ 
ticularly clear for the normal-stress difference. While its 
evolution after startup displays overshoots for large twain 
the decays after cessation is always monotonic, dictated 
by the fact that axx — (Tyy > 0 always holds in the model 
for arbitrary time-dependent laminar channel flow. 

To demonstrate the relevance of twaii for the initial 
startup and cessation evolution, we show in Fig. the 
velocity and shear-stress transients presented in Fig. 
as functions of t/twaii- For a Newtonian fluid, all curves 
for different channel widths would collapse. The same 
holds for the instantaneous nonlinear Maxwell model, 
Eqs. (15), discussed in Ref. [51 These results are shown 
as dashed lines in Fig. |^for comparison. Curves for the 
full nonlinear Maxwell model do not collapse, since struc¬ 
tural relaxation introduces a separate time scale that is 
independent on twain 


As shown in the right panels of Fig. the decay of 
both the velocity and the shear stress towards zero af¬ 
ter removal of the pressure gradient is oscillatory. The 
oscillations in the two quantities are shifted in phase, fol¬ 
lowing an initially faster decay of the velocity. This is a 
consequence of the history integral appearing in Eq. (|^ : 
as the velocity and with it the velocity gradients decay, 
the integral determining the stresses is still dominated 
by past contributions. Eventually, the velocity decays 
to zero, while stresses are still present. To relax these 
stresses, the fluid continues to flow, but in a direction 
opposite to the initial steady state (i.e., with negative 



t/t-vaH V^wall 

FIG. 3. Startup (left) and stopping (right) flow for the non¬ 
linear Maxwell model for different channel diameters, as in 
Fig. The velocity u{t) (top) and shear stress o'xy{i) (bot¬ 
tom) are shown as functions of t/twaii. Dashed lines represent 
the instantaneous Maxwell model, cf. text. 
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FIG. 4. Evolution of velocity (left axis) and shear stress (right 
axis) as a function of f/fwaii after cessation of pressure driven 
channel flow, for a channel corresponding to fwaii =???, at 
position y/H = 0.51. Solid lines are results for the nonlin¬ 
ear Maxwell model, dashed lines represent the instantaneous 
model. 


velocity since it is now driven by the remaining internal 
stresses rather than the external pressure gradient). This 
counter-flow causes the stresses to relax and eventually 
become negative, such that the velocity starts to increase 
towards positive values again. 

Figure [^highlights this phase-shifted oscillatory decay 
of the velocity and the shear stress for a single channel 
width. It is instructive to compare the observed decay 
pattern to that predicted by the instantaneous model 
(dashed lines in Fig. j^. Here, no oscillations are ob¬ 
served. Up to the first zero crossing of the velocity seen 
for the integral model, the decay of the velocity and the 
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FIG. 5. Startup velocity as in Fig. for channel widths 
corresponding to twaii/To = 0.1, 1, and 4, in the integral 
nonlinear Maxwell model with different structural relaxation 
times {8 = 10: dashed; 8 — 100: solid), evaluated at position 
y/H = 0.51. 


stress are very similar in the instantaneous model. There, 
however, the stresses quickly decay to zero once the ve¬ 
locity and its gradients vanish. The oscillatory behavior 
seen in the integral model is hence a true signature of 
viscoelasticity. 

For instantaneous models describing yield-stress fluids, 
such as Eqs. (151 in the limit r —>■ oo, a finite stopping 
time is observed in the decay of the velocitj®lllIl at some 
time tc, the stresses in the channel have all decayed to 
values below the yield stress, and hence the velocity has 
to obey u{t) = 0 exactly for all t > tc- Our model with 
T = IOtq does not have a true yield stress, but never¬ 
theless the instantaneous model shows the signature of 
this finite-time singularity. Increasing t, the kink visi¬ 
ble in the velocity decay around t/twaii = 1-2 becomes 
more pronounced. Realistic yield-stress fluids will typi¬ 
cally be viscoelastic, since the emergence of a yield stress 
is usually coupled to slow structural relaxation and its 
modification through the flow. In these fluids, the finite¬ 
time singularity tc does not mark the exact coming to 
rest of the flow, but rather sets a typical time scale for 
the oscillatory decay of velocities and stresses. 

We briefly discuss the influence of the structural re¬ 
laxation time on the transients. Figure compares the 
startup velocities for selected twaii (also shown in Fig. 
for r = IOtq and r = lOOro. These curves differ essen¬ 
tially only by the different steady-state values the tend to. 
For the smaller r, the steady-state velocities are higher, 
as the fluid has a lower viscosity in the low-shear region 
near the center of the channel. Still, for twaii ^ tq, an 
overshoot is seen that vanishes for twaii ^ ''"o- Note that 
in our definition of twaih the structural relaxation time r 
does not enter. This may appear surprising, since for vis¬ 
coelastic models, one might expect a strong dependence 
of the transient dynamics on the structural relaxation 
time. 



FIG. 6. Starting flow of the velocity in a narrow chan¬ 
nel (left, twaii/ro = 0.1) and an intermediate-width chan¬ 
nel (right, twaii/To = 1), evaluated at positions y/H = do 
and di /2 marked in Fig. [T| Solid lines represent the nonlin¬ 
ear Maxwell model including shear thinning, dashed lines the 
linear-viscoelastic UGM model, both for 8 = 10. Dash-dotted 
and dotted lines represent the instantaneous models with and 
without shear thinning (Eq. (151 and a Newtonian fluid, re¬ 
spectively) . 


To highlight the separate effects of (linear) viscoelastic¬ 
ity and shear thinning, we compare in Fig. the startup 
velocities for the integral nonlinear Maxwell model with 
those for the purely linear-viscoelastic UCM for. The lat¬ 
ter does not include shear thinning, so that the steady- 
state velocities are much lower, owing to the high vis¬ 
cosity set by the structural relaxation time r. The tran¬ 
sients observed in the UCM display overshoots (and, in 
fact, oscillations) for all channel widths, as highlighted 
for both twaii ^ To and twaii = tq in the figure. The 
inclusion of shear thinning into this viscoelastic model 
causes the oscillations to disappear. Realistic viscoelas¬ 
tic fluids will most likely also show shear thinning, since 
the appearance of slow structural relaxation makes the 
system prone to exhibit nonlinear-response phenomena. 
The qualitative behavior of the transient flow dynamics 
should hence be closer to our nonlinear model than to 
the UCM. 


C. Profile Evolution 

So far, we have discussed the time evolution of veloci¬ 
ties and stresses at selected positions across the channel. 
For a Newtonian fluid or the instantaneous modeP, this 
contains the essential information, since the shape of the 
cross-channel profiles does not change qualitatively dur¬ 
ing startup or cessation of the flow. 

The profile evolution in particular during cessatino of 
the flow is more complex for the present integral nonlin¬ 
ear Maxwell model. As shown in Fig. the velocity first 
relaxes to zero at slightly different times, depending on 
the cross-channel position. Interestingly, the plug in the 
center of the channel does not come to rest as a plug, but 
rather the velocity around y/H — 0 decreases faster than 
the nearby velocities. For the UCM without shear thin¬ 
ning, but including the Newtonian high-shear viscosity 
? 7 oo, this is not observed, and rather the center-channel 
velocity is slower to decay. 


























FIG. 7. Stopping flow of the velocity for the narrow (top, 
twaii/To = 0.1) and intermediate (bottom, twaii/ro = 1) chan¬ 
nel, for the integral nonlinear Maxwell model. The profiles 
(right) are plotted for different times in intervals indicated by 
a horizontal line of the same color at the top in the left panels. 
Each line of the same color is separated by At, which is dou¬ 
bled with each new color starting with = O.ltwaii 

and O.OStwaii, respectively. Profiles plotted with bold lines are 
taken at times marked by vertical lines in the left panels. 


The decay pattern of the velocities transmits to an 
intricate decay pattern also for the normal-stress differ¬ 
ence. Recall that we are considering laminar, translation- 
invariant, incompressible flow. The normal stresses are 
hence determined by the velocity gradients, but do not 
couple back to the flow evolution. 

Figure shows the evolution of the normal-stress- 
difference profiles in startup flow, for the two channel 
widths discussed also in Fig. These profiles evolve 
towards the characteristic steady-state profile which is 
quadratic in the center of the channel, crossing over to a 
constant near the walls. In the wider channel (fwaii = "T)), 
the increase in normal-stress difference close to the wall 
is more pronounced, and causes transient profiles that 
display a minimum at y/H = 0, and a maximum be¬ 
tween y/H = 0.5 and 1. In the cuts at constant y/H 
discussed above (see also left panel of the figure), this 
manifests itself as an overshoot in the time evolution 
which is not present for the narrower channel. Hence 
this overshoot is dominated by effects close to the chan¬ 
nel wall, where high shear rates transiently cause large 
normal-stress differences. Note that small compressibil¬ 
ity effects and normal-stress-induced particle migration 
might change this behavior. 

As mentioned above, in the cessation flow, normal- 
stress differences decay monotonically since they have to 
remain positive at all times. The corresponding profiles 
are shown in Fig. Again, the normal-stress difference 
shows the fastes transient evolution near the walls. Since 
in the center of the channel, the normal-stress difference 
is close to zero even in steady-state, this results in pro¬ 



0 20 40 60 80 100 120 140 -1 - 0.5 0 0.5 1 



0 2 4 6 8 10 12 14 -1 - 0.5 0 0.5 1 

V^waii y/H 

FIG. 8. Startup flow of the normal stress difference for the 
narrow (top) and intermediate (bottom) channel. The profiles 
(right) are plotted for different times indicated by a horizon¬ 
tal line of the same color at the top (left). Each line of the 
same color is separated by At, which is doubled with each 
new color starting with = O.ltwaii and Itwaii, re¬ 

spectively. Profiles plotted with bold lines are taken at times 
marked by vertical lines, the dashed line is the steady state 
profile. 

files that again have a minimum around y/H = 0 and 
a maximum at intermediate y/H. Different from what 
is seen in Fig. for the startup flow, in cessation, no 
qualitative change is observed between the narrow and a 
wider channel. 


V. CONCLUSION 

We have developed a hybrid-lattice-Boltzmann sim¬ 
ulation scheme that allows to simulate the flow of 
non-Newtonian, glass-forming fluids incorporating flow- 
history effects that arise from slow structural relaxation. 
The scheme builds upon an extension of the standard 
LB scheme to non-Newtonian constitutive equations pre¬ 
sented in Ref.l^. Here we extend it to include an integral- 
equation solver adapted to constitutive equations of the 
form, Eq. generically expected in nonlinear glassy 
rheology. The scheme is particularly adapted to deal with 
flows that include long-lived memory effects. 

The hybrid-LB algorithm was used to study the com¬ 
bined effects of viscoelasticity and shear thinning in 
pressure-driven planar channel flow of an incompressible 
fluid. To mimic features expected from microscopic the¬ 
ory, such as ITT-MCT^, we have employed a nonlinear 
generalized Maxwell model. The steady-state profiles of 
this model have been discussed earlier. As is typical for a 
fluid close to the glass transition, plug-like flow develops 
in the center of the channel, as a signature of the yield 
stress that arises at the glass transition. 

The transient evolution of velocities, shear stresses. 
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0 2 4 6 8 in 12 14 -1 -0.5 0 0.5 1 



0 1 2 3 4 5 -1 -0.5 0 0.5 1 

y/^ 

FIG. 9. Stopping flow of the normal stress difference for the 
narrow (top) and intermediate (bottom) channel. The profiles 
(right) are plotted for different times indicated by a horizontal 
line of the same color at the top (left). Each line of the same 
color is separated by At, which is doubled with each new color 
starting with — Q.ltwaii and O.OStwaii, respectively. 

Profiles plotted with bold lines are taken at times marked by 
vertical lines. 




FIG. 10. Startup (left) and stopping flow (right) for the UCM 
model for two channel diameters with twaii = O.Itq (squares) 
and twaii = To (circles). The velocity is measured at d = 0.51, 
lines are analytic results and symbols are obtained from LB 
simulations. 


equations with long-lived memory kernels. Such a com¬ 
bined MCT-LB scheme will be the subject of future work. 
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and normal-stress differences is rich in phenomenology. 
Since the time scale twaii characterizing hydrodynamic 
momentum transport across the channel is independent 
on the structural relaxation phenomena, the transient 
evolution of the flow differs in narrow channels from 
that in wide channels. For channel widths correspond¬ 
ing to twaii ^ To, overshoots appear in the startup veloc¬ 
ity. Note that twaii ~ tq corresponds to channel widths 
H = 0(mm), so that the effects should readily be ob¬ 
servable in experiment, e.g., on colloidal suspensions. 

The decay of the velocities and stresses after removal of 
the driving pressure gradient is oscillatory, reflecting the 
viscoelasticity that causes the stress evolution to lag be¬ 
hind that of the velocity. Transiently, the system comes 
to a rest, at a time given by the finite-time singular¬ 
ity discussed for yield-stress fluids with an instantaneous 
relation between stress and strain rate. As at this time, 
the stress has not fully decayed, it causes a backward mo¬ 
tion of the fluid, and hence slower oscillatory approach 
to rest. Such oscillations are already expected from lin¬ 
ear viscoelasticity, as exemplified by the upper-convected 
Maxwell model. However, generically, glass-forming flu¬ 
ids will exhibit both viscoelasticity and shear thinning. 

The nonlinear Maxwell model implemented here 
should account for many qualitative effects expected 
from more microscopic theories, such as ITT-MCT. Our 
hybrid-LB algorithm is readily adapted to include consti¬ 
tutive equations directly taken from ITT-MCT, although 
they are numerically much more demanding, since the ad- 
hoc exponential assumed in Eq. (§ is replaced by an ex¬ 
pression evolving density-correlation functions that need 
to be calculated from the solution of integro-differential 


Appendix A: Upper-convected Maxwell model 

The 2D Poiseuille flow of a UCM fluid has been solved 
analyticalljEH and provides a good test for our scheme to 
implement viscoelasticity in lattice Boltzmann simula¬ 
tions via an integral constitutive equation. Fig. [T0| shows 
the center velocity {d = y/H = 0.51) for two different 
channel diameters with twaii = O.Itq (red) and twaii = tq 
( blue) after applying (left) and removing (right) a sudden 
pressure difference. The channel diameters are the ones 
we find most interesting when discussing the nonlinear 
Maxwell model. Although only 100 nodes in transverse 
flow direction were used, the LB simulations reproduce 
the analytic results (black lines) extremely well. Devia¬ 
tions are larger in amplitude than in time and a higher 
precision can be easily attained by increasing the lattice 
size. The algorithm shows the same precision for the 
stopping flow as under startup and reproduces the sym¬ 
metry. 

To further test the capabilities of the LB scheme, we 
consider highly viscous UCM fluids in very wide chan¬ 
nels, e = {400, 2000} and 2H = L = {V|, V^jm. This 
way, the dimensionless retardation time^Sl S 2 = Goo(r -|- 
t‘o)t‘o/(pL^) = {0.05,0.01} <C 1, but the relaxation time 
Si = Goo(r-b ro)r/(pL 2 ) = {20.05,20.01} > 1. The 
density is p = lOOOkg/m^, the shear modulus Gao = IPa- 
This limit is interesting as the UCM fluid behaves similar 
to a soft elastic solid, but the large channel diameter al¬ 
lows the material to deform for a long time unperturbed 
by boundary effects. Fig [T^ shows the evolution of the 
center and half-center velocity in time. Please note, that 
the time is now given in units of twanTo/(T + tq) oc 
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+ To)/To i/«wau(T + Tq) / Tq 

FIG. 11. UCM. Velocity (solid lines) and shear stress (dashed 
lines, multiplied by 10) for 0 — 400 (left) and 6 = 2000 (right). 
The velocity and shear stress are scaled by their maximum 
steady state value. 

where jymax = Goo{,t + tq) is the long time limit of 
the viscosity under steady shear. Previously, the tran¬ 
sient dynamics was accelerated by shear and we identi¬ 
fied twaii oc T]^ as the characteristic transient time scale 
of the nonlinear Maxwell model. 

The initial response of the UCM model is almost purely 
elastic. Fig. 11 shows the evolution of the velocity (solid 
lines) and shear stress (dashed lines) in time at equally 
spaced positions in the channel. Red lines are in the chan¬ 
nel center, green lines at half-center. The applied pres¬ 
sure gradient exerts a constant body force on the fluid. 
As viscous damping takes place on a time scale much 
larger than the elastic response, the velocity increases 
linearly from each wall to develop a homogeneous shear 
field. Once the shear waves meet in the channel center, 
the fluid slows down again. The shear stress only starts 
to build up, when the velocity gradient is almost con¬ 
stant. It then increases linearly to satisfy the constant 
stress to strain relation of the Maxwell model. The driv¬ 
ing force is in turns used to increase either the kinetic 
energy or stress of the fluid. 

Viscous damping forces are small initially and only af¬ 
fect the dynamics on long time scales. As 9 is large, the 
memory of the initial state only decays slowly to finally 
give a flowing steady state much slower than the initial 
response after applying the pressure gradient, see Fig.[T3| 
The lattice Boltzmann scheme is able to track this long 
time evolution even for small lattice sizes. More precise 
results, especially for the first period, can be obtained 
when a larger lattice is used. We find the largest devi¬ 
ations from the analytic solution when there are sudden 
changes in the velocity due to the solid-like dynamics, as 
the LB algorithm always assumes a fluid. 
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